function [DSurplus,DSurplus_L,DSurplus_H,Oaxaca]=SurplusChangeWoman(GAMMA,STORE,STORE_,STORE1,STORE1_,Weight,natio,singl)
global dimS dimProv dimNum
% auxG=[GAMMA';1-GAMMA'];
% auxG=kron(auxG,ones(dimS,1));
% auxG=reshape(auxG,2*dimS*dimProv*dimNum,1);
auxG=[repmat(GAMMA'.*kron(Weight,ones(1,dimNum)),[dimS 1]);repmat((1-GAMMA').*kron(Weight,ones(1,dimNum)),[dimS 1])];
auxG=auxG(:);
auxSurp0=[STORE;STORE_];
% Tracking the identity of the woman across periods and markets.
e1=repmat((1:dimNum*dimProv)/1000,[dimS*2  1]);e1=permute(e1,[1 3 2]);
e2=repmat(kron([0.1;0.4],ones(dimS,1)),[1 1 dimNum*dimProv]);
e=e1+e2;
auxSurp0(:,8,:)=auxSurp0(:,8,:)+e;
auxSurp0=permute(auxSurp0,[1 3 2]);auxSurp0=reshape(auxSurp0,2*dimS*dimProv*dimNum,14);
if singl==-1
    if natio>0
    ind=(auxSurp0(:,1)==natio);
    elseif natio<0
      ind=(auxSurp0(:,1)~=abs(natio));  
    elseif natio==0
      ind= true(2*dimS*dimProv*dimNum,1);  
    end
else
    if natio>0
    ind=(auxSurp0(:,1)==natio)&(auxSurp0(:,7)==singl);
    elseif natio<0
    ind=(auxSurp0(:,1)~=1)&(auxSurp0(:,1)~=abs(natio))&(auxSurp0(:,7)==singl);
    elseif natio==0
    ind=logical(ones(2*dimS*dimProv*dimNum,1));  
    end
end
auxSurp1=[STORE1;STORE1_];
auxSurp1(:,8,:)=auxSurp1(:,8,:)+e;
auxSurp1=permute(auxSurp1,[1 3 2]);auxSurp1=reshape(auxSurp1,2*dimS*dimProv*dimNum,14);

Identity0=auxSurp0(ind,8);
ss0=zeros(length(Identity0),1);
ss1=zeros(length(Identity0),1);
ss0_L=zeros(length(Identity0),1);
ss0_H=zeros(length(Identity0),1);
ss1_L=zeros(length(Identity0),1);
ss1_H=zeros(length(Identity0),1);
Drho=zeros(length(Identity0),1);
DS1=zeros(length(Identity0),1);
DS2=zeros(length(Identity0),1);
DS3=zeros(length(Identity0),1);
DS4=zeros(length(Identity0),1);
ww=zeros(length(Identity0),1);
parfor i=1:length(Identity0)
    id0=find(auxSurp0(:,8)==Identity0(i));
    id1=find(auxSurp1(:,8)==Identity0(i));
    ss0_L(i)=auxSurp0(id0,10);
    ss0_H(i)=(1-auxSurp0(id0,7)).*(auxSurp0(id0,11)-auxSurp0(id0,9))+auxSurp0(id0,7).*auxSurp0(id0,10);
    ss0(i)=(1-auxSurp0(id0,7)).*auxSurp0(id0,11).*auxSurp0(id0,10)./(auxSurp0(id0,9)+auxSurp0(id0,10))+auxSurp0(id0,7).*auxSurp0(id0,10);
    ss1_L(i)=auxSurp1(id1,10);
    ss1_H(i)=(1-auxSurp1(id1,7)).*(auxSurp1(id1,11)-auxSurp1(id1,9))+auxSurp1(id1,7).*auxSurp1(id1,10);
    ss1(i)=(1-auxSurp1(id1,7)).*auxSurp1(id1,11).*auxSurp1(id1,10)./(auxSurp1(id1,9)+auxSurp1(id1,10))+auxSurp1(id1,7).*auxSurp1(id1,10);
    Drho(i)=auxSurp1(id1,10)./(auxSurp1(id1,9)+auxSurp1(id1,10))-auxSurp0(id0,10)./(auxSurp0(id0,9)+auxSurp0(id0,10));
    DS1(i,1)=(1-auxSurp0(id0,7)).*(1-auxSurp1(id1,7)).*auxSurp1(id1,11)*Drho(i); % Stay married, change in weight
    DS2(i,1)=(1-auxSurp0(id0,7)).*(1-auxSurp1(id1,7))*auxSurp0(id0,10)./(auxSurp0(id0,9)+auxSurp0(id0,10))*(auxSurp1(id1,11)-auxSurp0(id0,11)); % D surplus
    DS3(i,1)=auxSurp0(id0,7).*(1-auxSurp1(id1,7))*(ss1(i)-auxSurp0(id0,10)); % Single to married
    DS4(i,1)=(1-auxSurp0(id0,7)).*auxSurp1(id1,7)*(auxSurp1(id1,10)-ss0(i)); % Married to single
    ww(i)=auxG(i);
end
ww=ww/sum(ww);
S0=[ss0_L ss0_H ss0_L ss0_H];
S1=[ss1_L ss1_L ss1_H ss1_H];
DSurplus=S1-S0;
[DSurplusMin,i2]=min(DSurplus,[],2);
i3=sub2ind(size(DSurplus),(1:length(Identity0))',i2);
[DSurplusMax,i2]=max(DSurplus,[],2);
i4=sub2ind(size(DSurplus),(1:length(Identity0))',i2);

DSurplus_L=sum(DSurplusMin./S0(i3).*ww);
DSurplus_H=sum(DSurplusMax./S0(i4).*ww);
DSurplus=sum((ss1-ss0)./ss0.*ww);
% Oaxaca decomposition
Oaxaca=sum(ww.*[DS1 DS2 DS3 DS4]./ss0);
%Oaxaca=Oaxaca/sum(Oaxaca)*100;